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ABSTRACT 





Context. This paper is concerned with the image reconstruction problem when the mea- 
sured data are solar hard X-ray modulation profiles obtained from the Reuven Ramaty 
High Energy Solar Spectroscopic Imager (RHESSI) instrument. 

Aims. Our goal is to demonstrate that a statistical iterative method classically applied 
to the image deconvolution problem is very effective when utilized for the analysis of 
count modulation profiles in solar hard X-ray imaging based on Rotating Modulation 
Collimators. 

Methods. The algorithm described in this paper solves the maximum likelihood prob- 
lem iteratively and encoding a positivity constraint into the iterative optimization 
scheme. The result is therefore a classical Expectation Maximization method this time 
applied not to an image deconvolution problem but to image reconstruction from count 
modulation profiles. The technical reason that makes our implementation particularly 
effective in this application is the use of a very reliable stopping rule which is able to 
regularize the solution providing, at the same time, a very satisfactory Cash-statistic 
(C- statistic). 

Results. The method is applied to both reproduce synthetic flaring configurations and 
reconstruct images from experimental data corresponding to three real events. In this 
second case, the performance of Expectation Maximization, when compared to Pixon 
image reconstruction, shows a comparable accuracy and a notably reduced computa- 
tional burden; when compared to CLEAN, shows a better fidelity with respect to the 
measurements with a comparable computational effectiveness. 

Conclusions. If optimally stopped, Expectation Maximization represents a very reli- 
able method for image reconstruction in the RHESSI context when count modulation 
profiles are used as input data. 



Key words. Methods: Statistical image reconstruction — Methods: Expectation 
Maximization — Methods: regularization — Sun: flares 



2 Benvenuto et al.: Expectation Maximization for hard X-ray modulation profiles 

1. Introduction 



Expectation Maximization (EM) ([Dempster et al 19771) is an iterative algorithm that ad- 
dresses the maximum likelihood problem when 1) the relation between the unknown pa- 
rameter (or set of parameters) and the measured data is linear; 2) the data are drawn 
from a Poisson distribution; 3) the unknown parameter satisfies the positivity constraint. 
EM represe nts a genera lization of the image restoration method introduced by Lucy and 
Richardson (jLucv 19741) in the case of deconvolution problems typical of focused astron- 
omy and is successfully applied to several reconstruction problems in optics, microscopy 
and medical imaging. The present paper applies EM for the first time to the hard X-ray 
count modulation profiles measured by the nine Rotating Modulation Collimators (RMCs) 
mounted on th e Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) 
(jLin et al 20021) . More specifically, we show here that, when combined with an optimal 
stopping rule for the iterative process, EM provides reliable reconstructions with notable 
computational effectiveness. 

The RHESSI imaging concept translates to a s olar context the rota tional modulation 
synthesis first introduced for non-solar observations (jHurford et al 2 002). In RHESSI, a set 
of nine rotating collimators, characterized by nine pairs of grids with nine different pitches, 
time-modulates the incoming photon flux before it is detected by the nine corresponding 
Ge crystals. The resulting signal is a set of nine time series representing the count evolution 
provided by each collimator-detector system at different time bins. Therefore, these count 
modulation profiles represent the temporal or rotation angle/phase variation of the count 
rates for each grid. 

Since the transformation from the flux distribution on the image plane to the set of 
count modulation profiles is linear, the RHESSI image reconstruction problem is the linear 
inverse problem of describing the flux distribution from the count modulation profiles. EM 
describes the observed data and the unknown as realizations of stochastic quantities and 
searches for the flux distribution that maximizes the probability of the observation under 
the constraint that the pixel content must be positive. In fact, when the noise distribution 
on the data is Poisson this constrained maximum likelihood problem can be transformed 
into a fixed point problem whose solution is obtained iteratively, by means of a successive 
approximation scheme. From a theoretical viewpoint, the convergence properties of this 
algorithm when the number of iterations grows are not completely known. However, it 
is always observed in applications that stopping the procedure at some optimal iteration 
regularizes the reconstruction, thus preventing the occurrence of over-resolving effects or 
small-wavelength artifacts. In this paper, this optimized stopping rule is determined by 
utilizing the concept of constrained residual and by imposing that the empirical expectation 
value of this stochastic variable coincides with its theoretical expectation value. 

The plan of the paper is as follows. In Section 2 we describe the RHESSI imaging 
concept in more detail and introduce the linear transformation modeling the data formation 
process. Section 3 contains the formulation of the iterative algorithm together with its 
stopping rule. Section 4 validates EM in the case of several synthetic modulation profiles 
simulated from plausible configurations of the flux distribution. Finally, in Section 5 we 
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consider the application of EM to three sets of observed RHESSI data. Our conclusions 
are in Section 6. 

2. RHESSI count modulation profiles 

The RHESSI imaging hardware is made of nine sub-collimators, each one consisting of 
a pair of separated grids in front of a hard X-ray / gamma-ray detector. In each sub- 
collimator the two grids are identical, parallel and characterized by a planar array of 
equally-spaced, X-ray-opaque slats separated by transparent slits. The grid pitches of the 
different collimators are arranged according to a geometric progression with factor V3 
with detector 1 providing the maximum resolution power and minimum signal-to-noise 
ratio and detector 9 providing the minimum resolution power and maximum signal-to- 
noise ratio. RHESSI rotates around its own axis with a period of around 4 s and this 
rotation, combined with the presence of the grids, induces the modulation of the count 
rates. Therefore in this framework, there is no detector plane containing physical pixels 
(like in focused imaging) and here pixels are just a mathematical idealization in the image 
reconstruction process. If we describe the brightness distribution by means of the vector / 
of dimension iV 2 x 1 (in this lexicographic ordering each vector index denotes one of the 
N 2 pixels), then the expected count modulation detected by sub-collimator I is given by 

9 {l) = H®f, (1) 

where the vector has dimension pW x 1, pW is the number of time bins discretizing 
the time evolution of the modulation for detector I and H^\ with dimension pW x TV 2 , 
is the matrix modeling the transformation from the image to the measurement space. The 
entries of H^ can be interpreted in a probabilistic way as 

= AP®AU (2) 

where A measures the detector area and vf^ is the probability that a photon originating 
in pixel m will be counted in the z-th time bin of detector I during the time interval Ai$. 
If the image analysis performed involves data collected by M of the nine sub-collimators, 
then the overall signal formation process is described by 

.9 = Hf, (3) 

where g has dimension L x 1, L = J^fti an d contains all the modulation profiles while 
H has dimension L x N 2 and represents the image operator mimicking the action of the 
telescope. The EM algorithm together with its optimal stopping rule provides an estimate 
of / by means of an iterative regularized inversion of H. 

3. The EM algorithm 

Expectation Maximization is a statistical algorithm maximizing the probability that the 
data vector is a realization corresponding to a Poisson random vector g. In fact, in this 
case, the likelihood, i.e. the probability to observe g from the model Hf, can be written as 

i=l 9l - 
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where the ratio should be intended point-wise, element by element. We observe that 
maximizing this probability corresponds to minimizing the Cash statistic (C-statistic) 
jcash 197<Jl 

2 L 

c stat (g, V = 9i l0 § j^jy. + ( H f)i ~ 9i- (5) 

The likelihood maximizer is constrained to the set of positive solutions, i.e. the algorithm 
solves the constrained optimization problem 

argmaxP( 5 |/) = argminCW^, /), (6) 

which can be transformed into a fixed-point problem solved by means of the successive 
approximation scheme 

H T (^) 

Zfc+l = fk JjTi ' (^) 

with a positive (constant) initialization and where 1 denotes the vector made of all unit 
entries. Since H is ill conditioned, this iterative algorithm should be regularized by applying 
some stopping rule. To this aim we observe that the asymptotical behavior of equations 
([7]) is such that either fk — > or 

fe = — w\ — 

converges to 1. This implies that, asymptotically, 

tends to zero and therefore a reasonable stopping rule for EM in the RHESSI case is 

z k = E(z k ), (10) 
where E(zk) denotes the expectation value of Zf.- 



4. Numerical validation 

In order to assess the reliability of EM we setup a validation test based on the following 
process: 

1. Five different configurations of the flaring region were invented (see Figure [T]), first 
row, characterized by very different topographical and physical properties (e.g., size, 
position, number and distance of disconnected components, relative intensity of the 
components). Specifically, the original configurations are: a line source with constant 
density along the line (case A); a line source with intensity varying along the line, i.e. 
four compact sources and a weak one, all sources being aligned (case B); two Gaussian 
sources with flux ratio equal to 1 (case C); two Gaussian sources with flux ratio equal 
to 5 (case D); two Gaussian sources with flux ratio equal to 10 (case E); 

2. For each flaring configuration, three different synthetic sets of count modulation profiles 
were realized, characterized by three different levels of statistics (low, medium, high). 
Operationally, matrix H was applied to the simulated map, the resulting count expected 
values at each time bin were scaled with three different values in order to simulate three 
different levels of statistics (an average of 1000 counts per detector for the low level, 
10000 for the medium level and 100000 for the high one); 
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3. EM was applied to each one of the resulting 15 data sets in order to reconstruct the 
images; 

4. A set of routines was applied, for the quantitative assessment of the algorithm per- 
formance. These routines compute specific physical and geometrical parameters in the 
images, that are particularly significant for the different configurations, and compare 
the values with the corresponding ground truth values in the simulation maps of Figure 
[1] first row. 

Figure \T\ rows 2 through 4, contains the reconstructions provided by EM for the three 
different levels of statistics and using the count modulation profiles provided by all nine 
RHESSI detectors. The assessment routines are applied to these maps and compute the 
following parameters: 

Case A (line source with constant intensity): 

— Al: orientation (ground truth: deg). 

— A2: number of reconstructed sources (ground truth: 1). While reconstructing a line 
source, most (if not all) imaging methods tend to break it up into a set of compact 
sources (this, particularly, occurs at low levels of statistics). The routine computes 
the number K of intersection knots between the reconstructed line profile and the 
straight line with the same orientation passing through the Full Width at Half Maximum 
(FWHM). The number of reconstructed sources is counted as K/2. 

— A3: length (ground truth: 20.2 arcsec). The routine computes the FWHM in the direc- 
tion of the orientation line. 

— A4: width (ground truth: 1.35 arcsec). The routine computes the FWHM in the direction 
orthogonal to the orientation line. 

Case B (line source with intensity varying along the line) 

— Bl: orientation (ground truth: deg). 

— B2: number of reconstructed sources at FWHM (ground truth: 4). Computed as in case 
A2. 

— B3: length (ground truth: 16.9 arcsec). Computed as in case A3. 

— B4: width (ground truth: 1.35 arcsec). Computed as in case A4. 

Case C (two sources with flux ratio 1) 

— CI: position of the first reconstructed source (ground truth: arcsec). The routine 
computes the distance between the peak of the first reconstructed source and the cor- 
responding simulated one. 

— C2: position of the second reconstructed source (ground truth: arcsec). The routine 
computes the distance between the peak of the second reconstructed source and the 
corresponding simulated one. 

— C3: separation of the reconstructed sources (ground truth: 20 arcsec). The routine 
computes the distance between the two peaks. 

— C4: orientation of the separation line (ground truth: deg). The routine computes the 
orientation of the line passing through the two reconstructed peaks. 
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— C5: flux ratio (ground truth: 1). For each simulated source, the routine computes the 
disk centered in correspondence with the peak and with radius such that 99% flux 
is within the disk. Then, in the reconstructed image, the routine computes the fluxes 
contained in the two disks and makes the ratio. 

For Case D (two sources with flux ratio 5) and Case E (two sources with flux ratio 10) 
the routines compute the same parameters as in Case C. 

We performed this analysis by comparing the parameters obtained by EM with 
the original simulation par ameters and with the ones obt ained by di fferent imagin g 



algorithms, namely: Pixon (IPuetter 199a iMetcalf et al 19961 ). CLEAN (jrlogbom 19741) 



Maximum Entropy (Bong et al 20061 ) . a forward- fit algorithm for visibilities and uv .smooth 



( Massone et al 20091 ) (for all algorithms we used all nine RHESSI detectors). In Table Q] 
we reported just the results provided by Pixon since, among the methods using modulation 
profiles as input, provided among the best results. The Pixon algorithm models the source 
as a superposition of circular sources (or pixons) of different sizes and parabolic profiles 
and looks for the one that best reproduces the measured modulations from the different 
detectors. This technique is genera lly considered as the most reliable one in providing the 



most accurate image photometry (jDennis and Pernack 20091) . but at the price of a very 



notable computational burden. In this experiment we configured the Pixon algorithm in 
Solar Software (SSW) according to an optimized procedure based on heuristic arguments. 

5. Application to real observations 

We applied EM to RHESSI observations recorded in correspondence with three real flaring 
events. We considered the flare of April 15, 2002 in the time interval 00:06:00 - 00:08:00 
UT and in the energy interval 12 - 14 keV; the flare of February 20, 2002 in the time 
interval 11:05:58 - 11:06:41 UT and in the energy interval 25 - 30 keV; the flare of July 
23, 2002 in the time interval 00:30:00 - 00:32:00 UT and in the energy interval 100 - 300 
keV. In all cases detectors 3 through 8 are used for the observations. The reasons for this 
choice are as follows: detector 2 has been characterized by malfunctions since the beginning 
of the RHESSI mission; detector 1 is characterized by a very small signal-to-noise ratio 
while the coarse information carried by detector 9 is not crucial for the reconstruction of 
these events. Figure [5] compares EM reconstructions with the ones provided by Pixon and 
CLEAN, while Table [5] contains the corresponding C-statistic for the six detectors employed 
in the analysis. In these ex periments CLEAN param eters have been chosen according to 



optimized heuristic recipes (jDennis and Pernack 20 09). According to these results, EM and 
Pixon are characterized by values of the C-statistic almost systematically close to 1 (and 
significantly smaller than the ones provided by CLEAN). This means the EM and Pixon can 
reproduce the data with comparable accuracy (significantly better than the one achieved 
by CLEAN), although EM reduces the computational time of up to a factor 4. Indeed, for 
April 15, 2002, the computational time is around 100 sec for EM and more than 400 sec for 
Pixon; for February 20, 2002, the computational time is around 50 sec for EM and almost 
220 sec for Pixon; for July 23, 2002, the computational time is around 125 sec for EM and 
more than 460 sec for Pixon. This same 4 to 1 scaling holds true independently of the kind 
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Fig. 1. Validation of EM with synthetic data. First row: the simulated configurations. Rows 
2 through 4: reconstructions provided by EM corresponding to count modulation profiles 
characterized by three different levels of statistics (second row: average of 100000 counts 
per detector; third row: average of 10000 counts per detector; fourth row: average of 1000 
counts per detector. 

of hardware used for the tests. We also observe that in the Pixon implementation we used 
for these experiments the computations of the time profiles and back projections are done 
more efficiently by using optimized combinations of the spatial variation as spatial sine and 
cosine patterns (annsec, annular-sector, implementation). Implementing EM according to 
this same representation will improve the computational gain provided by this algorithm 
of another factor 5. 

6. Conclusions 

This papers shows that Expectation Maximization can be effectively applied to recon- 
struct hard X-ray images of solar flares from the count modulation profiles recorded by 
the RHESSI mission. This method is an iterative likelihood maximizer with a positivity 
constraint, that explicitly exploits the fact that the noise affecting the measured data has 
a Poisson nature. We have utilized an optimal stopping rule that regularizes the algorithm, 
realizing an optimal trade-off between the C-statistic and the numerical stability of the 
reconstruction. We are aware that this test on C-statistic in not conclusive, since images 
affected by significant artifacts may reproduce the experimental data with great accuracy. 
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Table 1. Assessment of EM performances in the case of synthetic data; 10 5 , 10 4 and 10 3 
(counts per detector) indicate the three different levels of statistics considered in the test. 
The comparison is made with the parameters characterizing the simulated configurations 
(ground truth, last column). Case A concerns the line source with constant intensity (Figure 
[TJ first row, first panel): orientation (Al, deg); number or reconstructed sources (A2); 
length (A3, arcsec); width (A4, arcsec). Case B concerns the line source with intensity 
varying along the line (Figure [TJ first row, second panel): orientation (Bl, deg); number of 
reconstructed sources (B2); length (B3, arcsec); width (B4, arcsec). Case C concerns two 
sources with flux ratio 1 (Figure [TJ first row, third panel): position of the first reconstructed 
source (CI, arcsec); position of the second reconstructed source (C2, arcsec); separation 
(C3, arcsec); orientation of the separation line (C4, deg); flux ratio (C5). Case D (Figure 
[TJ first row, fourth panel) and Case E (Figure [TJ first row, fifth panel) as for Case C but 
the flux ratio 5 and 10, respectively. 



However low C-statistic values coupled with the positivity constraint can be considered as 
a diagnostic of reliable reconstructions. 

We have validated EM against synthetic count modulation profiles corresponding to 
challenging simulated configurations and characterized by three levels of statistics. Then 
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!) 







Fig. 2. Performance of EM in the case of real data observed by RHESSI. First row: EM 
reconstructions; second row: Pixon reconstructions; third row: CLEAN reconstructions. 
First Column: 15 April 2002 event; second column: 20 February 2002 event; third column: 
July 23 2002 event. 



we have applied the method against the RHESSI observations of three flaring events and 
compared the reconstructions with the ones provided by Pixon and CLEAN. The results of 
these experiments show that EM combines a reconstruction fidelity (in terms of C-statistic) 
comparable with the one provided by Pixon (which, however, is much more demanding 
from a computational viewpoint) with a computational efficiency comparable with the one 
offered by CLEAN (which, however, predicts the count modulation profiles by means of a 
significantly worse C-statistic). Our next step, which is currently under construction, will 
be to generalize this approach to the reconstruction of electron flux maps of the flaring 
region. Ele ctron flux maps o f solar flares can be already generated by hard X-ray count 
visibilities (jPiana et al 20071 ). We are currently working at an EM-based approach to the 
reconstruction of electron images, where the input data are the count modulation profiles 
and the imaging matrix to invert accounts for both the effects of the brcmsstrahlung cross- 
section and the Detector Response Matrix mimicking the projection from the photon to 
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Table 2. Performance of EM in the case of real data observed by RHESSI: comparisons 
of C-statistic provided by EM, Pixon and CLEAN. 

the count domain. The advantage of this approach with the respect to the visibility-based 
one should be that EM provides an analysis framework that is closer to the data, as we 
can model with a greater accuracy all of the detector effects. 
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